module matrix

contains

    subroutine matrix_inverse(n, A, B)
        integer, intent(in) :: n
        real, intent(inout) :: A(n,n)
        real, intent(out) :: B(n,n)

        integer :: i, j, k
        real :: temp_array(n)
        real :: factor
        logical :: flag

        do i = 1, n
            B(i,i) = 1.0
        end do

        do k = 1, n-1
            if(A(k,k) == 0) then
                flag = .false.
                do j = k+1, n
                    if(A(j,k) /= 0) then
                        temp_array = A(k,:)
                        A(k,:) = A(j,:)
                        A(j,:) = temp_array
                        temp_array = B(k,:)
                        B(k,:) = B(j,:)
                        B(j,:) = temp_array
                        flag = .true.
                        exit
                    end if
                end do
                if(flag == .false.) then
                    ! handle error
                    write(*,*) "matrix_inverse: error: singular matrix"
                    stop
                end if
            end if
            do i = k+1, n
                factor = A(i,k)/A(k,k)
                A(i,:) = A(i,:)-factor*A(k,:)
                B(i,:) = B(i,:)-factor*B(k,:)
            end do
            B(k,:) = B(k,:)/A(k,k)
            A(k,:) = A(k,:)/A(k,k)
        end do

        do k = n, 2, -1
            do i = k-1, 1, -1
                factor = A(i,k)
                A(i,:) = A(i,:)-factor*A(k,:)
                B(i,:) = B(i,:)-factor*B(k,:)
            end do
        end do

        return
    end subroutine matrix_inverse

end module matrix
